% ferry correlations Spearman or Pearson 
corrMethod='Pearson';%'Spearman';%

% --- data (example from ferry slider 50)----
Rdelays=delays_Ferry_Slider_50';
Rscores=scores_Ferry_Slider_50';
Rsubjects=subjects_Ferry_Slider_50';


% remove training set
 Rdelays_temp=Rdelays;
 Rscores_temp=Rscores;
 Rsubjects_temp=Rsubjects;
 Rdelays=Rdelays_temp(Rdelays_temp~=20 & Rdelays_temp~=4);
 Rscores=Rscores_temp(Rdelays_temp~=20 & Rdelays_temp~=4);
 Rsubjects=Rsubjects_temp(Rdelays_temp~=20 & Rdelays_temp~=4); 

 
% ---  Within subjects correlations and sem ---
people_id=unique(Rsubjects); % list of subjects ID
people_num=length(people_id); % length of subjects ID


individual_corrs=zeros(1, people_num); % this is an array to store within subject Spearman R
for i=1:people_num
    x=Rscores(Rsubjects==people_id(i));
    y=Rdelays(Rsubjects==people_id(i));
    individual_corrs(1,i)=corr(x',y','type',corrMethod).^2;
end
individual_corrs=individual_corrs(~isnan(individual_corrs));
meanIndividualCorrs=mean(individual_corrs);
individual_se=std(individual_corrs)./sqrt(people_num);
% Results are: 
% meanIndividualCorrs
% individual_se


% --- Pooled correlations with bootstrap error bars ---
PoolCorr=corr(Rdelays', Rscores','type',corrMethod).^2;
% results: PoolCorr

% Bootstrpa by resampling subjects with replacement 1000 times
PoolCorrBoot=zeros(1000,1);
for boot=1:1000
    tmpScore=[];
    tmpDelay=[];
    for i=1:people_num
        rnd=floor(rand()*people_num)+1; % find a random subject
        randSubj=people_id(rnd);
        tmpDelay=[tmpDelay, Rdelays(Rsubjects==randSubj)];
        tmpScore=[tmpScore, Rscores(Rsubjects==randSubj)];
    end
    PoolCorrBoot(boot)=corr(tmpDelay', tmpScore','type',corrMethod).^2;
end
PoolCorrBoot=sort(PoolCorrBoot);
% Results:
% low CI: PoolCorrBoot(50)
% high CI: PoolCorrBoot(950)


% --- Averaged correlations with boot strap error bars ---
%  bin the pooled data and average
pooledbin=zeros(20,1);
index=1;
for i=0:2:39
    pooledbin(index)=mean(Rscores(Rdelays>=i & Rdelays<i+2));
    index=index+1;
end
x=1:20;
AvgCorr=corr(x',pooledbin,'type',corrMethod).^2;
%results: AvgCorr

% Bootstrpa by resampling subjects with replacement 1000 times
AvgCorrBoot=zeros(1000,1);
for boot=1:1000
    tmpScore=[];
    tmpDelay=[];
    for i=1:people_num
        rnd=floor(rand()*people_num)+1; % find a random subject
        randSubj=people_id(rnd);
        tmpDelay=[tmpDelay, Rdelays(Rsubjects==randSubj)];
        tmpScore=[tmpScore, Rscores(Rsubjects==randSubj)];
    end
    index=1;
    for i=0:2:39
        pooledbin(index)=mean(tmpScore(tmpDelay>=i & tmpDelay<i+2));
        index=index+1;
    end
    x=1:20;
    AvgCorrBoot(boot)=corr(x',pooledbin,'type',corrMethod).^2;
end
AvgCorrBoot=sort(SpearmanAvgCorrBoot);
% results 
% AvgCorrBoot(50)
% AvgCorrBoot(950)



% --- AUC (area under curve) with bootstrap error bars --- 
AverageSizeCorrs=zeros(34,100);
for iter=1:100
    for j=2:35 
        pooledbin=zeros(20,1);
        index=1;
        for i=0:2:39
            tmpS=Rscores(Rdelays>=i & Rdelays<i+2);
            % ensure unique subjects pool:
            tmpSubj=Rsubjects(Rdelays>=i & Rdelays<i+2);
            c=[tmpSubj',tmpS'];
            c=sortrows(c,1);
            c(1,3)=1;
            c(2:end,3)=diff(c(:,1));
            tmpSubj=tmpSubj(c(:,3)>0); % these are all unique subjects
            % randomize tmpS order, and choose the first js several times like
            % boot strap, so that we get say 10 such estimates for each number...
            rnd=rand(length(tmpS),1);
            [~,Bsort]=sort(rnd); %Get the order of B
            tmpS=tmpS(Bsort);
            pooledbin(index)=mean(tmpS(1:min(j,length(tmpS))));
            index=index+1;
        end
        x=1:20;
        AverageSizeCorrs(j-1,iter)=corr(x',pooledbin,'type',corrMethod).^2; 
    end 
end
AUCgraph=mean(AverageSizeCorrs,2);
AUC=mean(mean(AverageSizeCorrs));
% results: AUC



% Compute error bars for AUC using resampling with replacement of subjects bootstrap
aucFerryboot=zeros(1000,1);
for boot=1:1000
    tmpScore=[];
    tmpDelay=[];
    tmpSubject=[];
    for i=1:people_num
        rnd=floor(rand()*people_num)+1;
        randSubj=people_id(rnd);
        tmpDelay=[tmpDelay, Rdelays(Rsubjects==randSubj)];
        tmpScore=[tmpScore, Rscores(Rsubjects==randSubj)];
        tmpSubject=[tmpSubject, Rsubjects(Rsubjects==randSubj)];
    end
    
  % compute AUC 
    for iter=1:100 % 100 iterations of sampling 
        for j=2:35 
            pooledbin=zeros(20,1);
            index=1;
            for i=0:2:39
                tmpS=tmpScore(tmpDelay>=i & tmpDelay<i+2);
                % ensure unique subjects pool:
                tmpSubj=tmpSubject(tmpDelay>=i & tmpDelay<i+2);
                c=[tmpSubj',tmpS'];
                c=sortrows(c,1);
                c(1,3)=1;
                c(2:end,3)=diff(c(:,1));
                tmpSubj=tmpSubj(c(:,3)>0); % these are all unique subjects
                % randomize tmpS order, and choose the first js several times like
                % boot strap, so that we get say 10 such estimates for each number...
                rnd=rand(length(tmpS),1);
                [~,Bsort]=sort(rnd); %Get the order of B
                tmpS=tmpS(Bsort);
                pooledbin(index)=mean(tmpS(1:min(j,length(tmpS))));
                index=index+1;
            end
            x=1:20;
            AverageSizeCorrs(j-1,iter)=corr(x',pooledbin,'type',corrMethod).^2; 
        end   
    end
    aucFerryboot(boot)=mean(mean(AverageSizeCorrs));  
end
aucFerryboot=sort(aucFerryboot);
% results
% confLow=aucFerryboot(50)
% confHigh=aucFerryboot(950)